clear all;

cd ..; cd ..; cd ..;
load('code_main_model\6_implied_w_r\computation_results\data.mat','sigma','V_backward','nmkt','pdf_B_forward',...
    'nb','ntheta','nz','gridb','gridtheta','gridz','nb_pdf','ntheta_pdf','gridb_pdf','gridtheta_pdf');
cd 'code_figure_table\maintext\Figure 12';
t_ini=1;
t_end=11;
V_ini=squeeze(V_backward(t_ini,:,:,:,:));
V_new=squeeze(V_backward(t_end,:,:,:,:));

pdf_ini=squeeze(pdf_B_forward(t_ini,:,:,:,:));
pdf_new=squeeze(pdf_B_forward(t_end,:,:,:,:));

welfare(1:nmkt,1)=nan;
V0(1:nmkt,1)=0;
V1(1:nmkt,1)=0;
for i_mkt=1:nmkt
    for i_z=1:nz
        for i_b=1:nb_pdf
            for i_theta=1:ntheta_pdf
                V0(i_mkt)=V0(i_mkt)+lin_interpo2(nb,ntheta,gridb,gridtheta,squeeze(V_ini(i_mkt,i_z,:,:)),gridb_pdf(i_b),gridtheta_pdf(i_theta)) ...
                    .*pdf_ini(i_mkt,i_z,i_b,i_theta);
                V1(i_mkt)=V1(i_mkt)+lin_interpo2(nb,ntheta,gridb,gridtheta,squeeze(V_new(i_mkt,i_z,:,:)),gridb_pdf(i_b),gridtheta_pdf(i_theta)) ...
                    .*pdf_new(i_mkt,i_z,i_b,i_theta);
            end
        end
    end
    V0(i_mkt)=V0(i_mkt)./sum(pdf_ini(i_mkt,:,:,:),'all');
    V1(i_mkt)=V1(i_mkt)./sum(pdf_new(i_mkt,:,:,:),'all');
    welfare(i_mkt)=(V1(i_mkt)/V0(i_mkt))^(1/(1-sigma))-1;
end